MODULE fmm_boxed_multipole_ints

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fmm_init_boxed_multipole_ints,      &
             fmm_free_boxed_multipole_ints,      &
             fmm_build_boxed_multipoles,         &
             fmm_pack_boxed_mpoles,              &
             fmm_contract_multipole_potential

   REAL(REALK), ALLOCATABLE :: ECoefX(:,:,:,:) 
   REAL(REALK), ALLOCATABLE :: ECoefY(:,:,:,:)
   REAL(REALK), ALLOCATABLE :: ECoefZ(:,:,:,:)

   REAL(REALK), ALLOCATABLE :: MIntX(:,:,:)
   REAL(REALK), ALLOCATABLE :: MIntY(:,:,:)
   REAL(REALK), ALLOCATABLE :: MIntZ(:,:,:)

   REAL(REALK), ALLOCATABLE :: MpoleX(:,:,:,:)
   REAL(REALK), ALLOCATABLE :: MpoleY(:,:,:,:)
   REAL(REALK), ALLOCATABLE :: MpoleZ(:,:,:,:)

   REAL(REALK), ALLOCATABLE :: CarMpole(:,:,:)
   REAL(REALK), ALLOCATABLE :: SphMpole(:,:,:)

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_init_boxed_multipole_ints(basis,MaxMul)

      IMPLICIT NONE
      TYPE(fmm_basis), INTENT(IN) :: basis
      INTEGER(INTK),   INTENT(IN) :: MaxMul

      INTEGER(INTK) :: ndim, MaxAngl, MaxSgm2
      MaxAngl = basis%maxangl 
      MaxSgm2 = basis%maxsgm2 

      ALLOCATE(ECoefX(0:MaxAngl, 0:MaxAngl, 0:MaxAngl*2, MaxSgm2))
      ALLOCATE(ECoefY(0:MaxAngl, 0:MaxAngl, 0:MaxAngl*2, MaxSgm2))
      ALLOCATE(ECoefZ(0:MaxAngl, 0:MaxAngl, 0:MaxAngl*2, MaxSgm2))

      ndim = 2*MaxAngl + MaxMul
      ALLOCATE(MIntX(-1:ndim, -1:MaxMul+1, MaxSgm2))
      ALLOCATE(MIntY(-1:ndim, -1:MaxMul+1, MaxSgm2))
      ALLOCATE(MIntZ(-1:ndim, -1:MaxMul+1, MaxSgm2))
!      ndim = MaxMul+1
!      ALLOCATE(MIntX(-1:ndim, -1:MaxMul+1, MaxSgm2))
!      ALLOCATE(MIntY(-1:ndim, -1:MaxMul+1, MaxSgm2))
!      ALLOCATE(MIntZ(-1:ndim, -1:MaxMul+1, MaxSgm2))

      ALLOCATE(MpoleX(0:MaxAngl, 0:MaxAngl, 0:MaxMul, MaxSgm2))
      ALLOCATE(MpoleY(0:MaxAngl, 0:MaxAngl, 0:MaxMul, MaxSgm2))
      ALLOCATE(MpoleZ(0:MaxAngl, 0:MaxAngl, 0:MaxMul, MaxSgm2))

   END SUBROUTINE fmm_init_boxed_multipole_ints

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_free_boxed_multipole_ints

      IMPLICIT NONE
      DEALLOCATE(ECoefX)
      DEALLOCATE(ECoefY)
      DEALLOCATE(ECoefZ)

      DEALLOCATE(MIntX)
      DEALLOCATE(MIntY)
      DEALLOCATE(MIntZ)

      DEALLOCATE(MpoleX)
      DEALLOCATE(MpoleY)
      DEALLOCATE(MpoleZ)

      IF (ALLOCATED(CarMpole)) DEALLOCATE(CarMpole) 
      IF (ALLOCATED(SphMpole)) DEALLOCATE(SphMpole) 

   END SUBROUTINE fmm_free_boxed_multipole_ints

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_build_boxed_multipoles(scheme,basis,lmax,sh_pairs)

      USE fmm_car_to_sph,     ONLY: fmm_init_car_to_sph,      &
                                    fmm_free_car_to_sph
      USE fmm_integral_utils, ONLY: fmm_get_prim_batch,       &
                                    fmm_build_Ecoef1,         &
                                    fmm_build_Ecoef2
      USE fefmm_box_utils,    ONLY: fefmm_init_box_utils,     &
                                    fefmm_get_boxes,          &
                                    fefmm_box_centre,         &
                                    fefmm_sh_pair_in_box

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN) :: scheme
      TYPE(fmm_basis),    INTENT(IN) :: basis
      INTEGER(INTK),      INTENT(IN) :: lmax
      TYPE(fmm_sh_pairs), INTENT(IN) :: sh_pairs(:)

!fixme
      CHARACTER(LEN=10), PARAMETER :: Name = 'multipoles'

      TYPE(fmm_prim_batch) :: batch(basis%MaxSgm2)
      CHARACTER(LEN=255) :: FBuf
      INTEGER(INTK) :: Ish, Jsh, NPrim, IAnglA, IAnglB
      INTEGER(INTK) :: i, ij, n, nbox
      REAL(REALK)   :: grain, box(3)
      LOGICAL       :: empty_box

      grain = scheme%grain
      CALL fefmm_init_box_utils(scheme%grain,sh_pairs)
      CALL fefmm_get_boxes(nbox)
      IF (nbox < 3) RETURN
      print *, 'Building raw multipole moments...'

      CALL fmm_init_car_to_sph(lmax)

      FBuf = TRIM(Name)//".fmm0"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

      DO n = 1, nbox
         box = fefmm_box_centre(n)
!print *, 'box:',n, box
         empty_box = .TRUE.

         DO ij = 1, SIZE(sh_pairs)

            IF (.NOT.fefmm_sh_pair_in_box(sh_pairs(ij),box)) CYCLE

            Ish = sh_pairs(ij)%I
            Jsh = sh_pairs(ij)%J

            ! Prepare batch of primitives
            CALL fmm_get_prim_batch(basis,Ish,Jsh,batch,NPrim)

            ! Build Hermite expansion E-coefficients
            IAnglA = basis%KType(Ish)
            IAnglB = basis%KType(Jsh)
            IF (IAnglA >= IAnglB) THEN
               CALL fmm_build_Ecoef1(batch,NPrim,IAnglA,IAnglB,   &
                                     ECoefX,ECoefY,ECoefZ)
            ELSE
               CALL fmm_build_Ecoef2(batch,NPrim,IAnglA,IAnglB,   &
                                     ECoefX,ECoefY,ECoefZ)
            END IF

            ! Set multipole origin and integration range 
            SELECT CASE (scheme%job_type)
               CASE (GFC_FMM)
                  DO i = 1, NPrim
                     batch(i)%PC(:) = batch(i)%P(:) - box(:)
                     batch(i)%hi(:) = -1d99    ! P - hi
                     batch(i)%lo(:) =  1d99    ! P - lo
                  END DO
               CASE DEFAULT
                  DO i = 1, NPrim
                     batch(i)%PC(:) = batch(i)%P(:) - box(:)
                     batch(i)%hi(:) = batch(i)%P(:) - (box(:)+grain*half)
                     batch(i)%lo(:) = batch(i)%P(:) - (box(:)-grain*half)
                  END DO
            END SELECT
            ! Build truncated Hermite multipole moment integrals
            CALL fmm_build_boxed_Mints(batch,NPrim,IAnglA+IAnglB,lmax)
            !!!!CALL fmm_build_Mints(batch,NPrim,IAnglA+IAnglB,lmax)

            ! Build spherical multipole integrals for Cartesian GTOs
            ! assuming final Fock matrix in spherical GTOs
            ! (only difference is we use spherical AO contraction coefs)
            ! When using these multipoles need Cartesian form of
            ! density and J_matrices
            CALL fmm_build_boxed_Mpole_ints(IAnglA,IAnglB,NPrim,lmax)
            CALL fmm_build_boxed_SphMpole(basis,batch,Ish,Jsh,NPrim,lmax)
            CALL fmm_store_boxed_SphMpole(basis,batch,Ish,Jsh,lmax,empty_box)

         END DO

         ! Mark end of box contents
         WRITE(IO2) -1,-1, 0,0, 0d0
         ! Store box centre
         WRITE(IO2) empty_box, box(:)

      END DO

      ! Mark end of file
      WRITE(IO2) -2,-2, 0,0, 0d0
      CLOSE(UNIT=IO2, STATUS='KEEP')

      CALL fmm_free_car_to_sph

   END SUBROUTINE fmm_build_boxed_multipoles

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_build_boxed_Mints(batch,NPrim,MaxAng,MaxMul)

!     o Evaluate 1D multipole-moment integrals

      IMPLICIT NONE

      TYPE(fmm_prim_batch), INTENT(IN) :: batch(:)
      INTEGER(INTK),        INTENT(IN) :: NPrim, MaxAng, MaxMul

      INTEGER :: i, Ie, It, IJ

      REAL(REALK) :: PXCX(NPrim), PYCY(NPrim), PZCZ(NPrim)
      REAL(REALK) :: amx(3,NPrim), bmx(3,NPrim)
      REAL(REALK) :: ExpPHalf(NPrim), ExpntP(NPrim)

      REAL(REALK) :: TaIntX(0:MaxMul+MaxAng, NPrim)
      REAL(REALK) :: TaIntY(0:MaxMul+MaxAng, NPrim)
      REAL(REALK) :: TaIntZ(0:MaxMul+MaxAng, NPrim)
      REAL(REALK) :: TbIntX(0:MaxMul+MaxAng, NPrim)
      REAL(REALK) :: TbIntY(0:MaxMul+MaxAng, NPrim)
      REAL(REALK) :: TbIntZ(0:MaxMul+MaxAng, NPrim)

      PXCX(:) = batch(1:Nprim)%PC(1)
      PYCY(:) = batch(1:Nprim)%PC(2)
      PZCZ(:) = batch(1:Nprim)%PC(3)
      ExpPHalf(:) = batch(1:Nprim)%ExpPHalf
      ExpntP(:) = batch(1:Nprim)%ExpntP

      DO i = 1, NPrim
         amx(:,i) = batch(i)%lo(:)
         bmx(:,i) = batch(i)%hi(:)
      END DO

!     o Recurence relation
!       M(t,e,NPrim)
!         e: 0 <= Ie <= MaxMul
!         t: 0 <= It <= MaxMul + MaxAng - e

      MIntX(:,:,:) = Zero
      MIntY(:,:,:) = Zero
      MIntZ(:,:,:) = Zero
      TaIntX(:,:) = Zero
      TaIntY(:,:) = Zero
      TaIntZ(:,:) = Zero
      TbIntX(:,:) = Zero
      TbIntY(:,:) = Zero
      TbIntZ(:,:) = Zero

!     o Initialise all the M_t for e=0

!     t=0, e=0 

      TaIntX(0,:) = ERF(SQRT(ExpntP(:))*amx(1,:))
      TaIntY(0,:) = ERF(SQRT(ExpntP(:))*amx(2,:))
      TaIntZ(0,:) = ERF(SQRT(ExpntP(:))*amx(3,:))

      TbIntX(0,:) = ERF(SQRT(ExpntP(:))*bmx(1,:))
      TbIntY(0,:) = ERF(SQRT(ExpntP(:))*bmx(2,:))
      TbIntZ(0,:) = ERF(SQRT(ExpntP(:))*bmx(3,:))

      MIntX(0,0,1:NPrim) = Half*SQRT(Pi/ExpntP(:)) * (TaIntX(0,:)-TbIntX(0,:))
      MIntY(0,0,1:NPrim) = Half*SQRT(Pi/ExpntP(:)) * (TaIntY(0,:)-TbIntY(0,:))
      MIntZ(0,0,1:NPrim) = Half*SQRT(Pi/ExpntP(:)) * (TaIntZ(0,:)-TbIntZ(0,:))

!     t=1, e=0 

      TaIntX(1,:) = EXP( -ExpntP(:)*amx(1,:)*amx(1,:) )
      TaIntY(1,:) = EXP( -ExpntP(:)*amx(2,:)*amx(2,:) )
      TaIntZ(1,:) = EXP( -ExpntP(:)*amx(3,:)*amx(3,:) )

      TbIntX(1,:) = EXP( -ExpntP(:)*bmx(1,:)*bmx(1,:) )
      TbIntY(1,:) = EXP( -ExpntP(:)*bmx(2,:)*bmx(2,:) )
      TbIntZ(1,:) = EXP( -ExpntP(:)*bmx(3,:)*bmx(3,:) )

      MIntX(1,0,1:NPrim) = TaIntX(1,:) - TbIntX(1,:)
      MIntY(1,0,1:NPrim) = TaIntY(1,:) - TbIntY(1,:)
      MIntZ(1,0,1:NPrim) = TaIntZ(1,:) - TbIntZ(1,:)

      DO It = 2, MaxMul+MaxAng

         TaIntX(It,:) = -Two*ExpntP(:)*( amx(1,:)*TaIntX(It-1,:)     &
                                         + (It-2)*TaIntX(It-2,:) )
         TaIntY(It,:) = -Two*ExpntP(:)*( amx(2,:)*TaIntY(It-1,:)     &
                                         + (It-2)*TaIntY(It-2,:) )
         TaIntZ(It,:) = -Two*ExpntP(:)*( amx(3,:)*TaIntZ(It-1,:)     &
                                         + (It-2)*TaIntZ(It-2,:) )

         TbIntX(It,:) = -Two*ExpntP(:)*( bmx(1,:)*TbIntX(It-1,:)     &
                                         + (It-2)*TbIntX(It-2,:) )
         TbIntY(It,:) = -Two*ExpntP(:)*( bmx(2,:)*TbIntY(It-1,:)     &
                                         + (It-2)*TbIntY(It-2,:) )
         TbIntZ(It,:) = -Two*ExpntP(:)*( bmx(3,:)*TbIntZ(It-1,:)     &
                                         + (It-2)*TbIntZ(It-2,:) )

         MIntX(It,0,1:NPrim) = TaIntX(It,:) - TbIntX(It,:)
         MIntY(It,0,1:NPrim) = TaIntY(It,:) - TbIntY(It,:)
         MIntZ(It,0,1:NPrim) = TaIntZ(It,:) - TbIntZ(It,:)

      END DO

!     t=2..(MaxMul+MaxAng) , e=0 
!
!     o Recurrences for remaining {t,e} same as usual, with different t-range
!
      DO Ie = 1, MaxMul
         DO It = 0, MaxMul+MaxAng-Ie
!         DO It = 0, Ie
            DO IJ = 1, NPrim
               MIntX(It,Ie,IJ) = DBLE(It) * MIntX(It-1,Ie-1,IJ)
               MIntY(It,Ie,IJ) = DBLE(It) * MIntY(It-1,Ie-1,IJ)
               MIntZ(It,Ie,IJ) = DBLE(It) * MIntZ(It-1,Ie-1,IJ)
            END DO
            DO IJ = 1, NPrim
               MIntX(It,Ie,IJ) = MIntX(It,Ie,IJ) + PXCX(IJ) * MIntX(It,Ie-1,IJ)
               MIntY(It,Ie,IJ) = MIntY(It,Ie,IJ) + PYCY(IJ) * MIntY(It,Ie-1,IJ)
               MIntZ(It,Ie,IJ) = MIntZ(It,Ie,IJ) + PZCZ(IJ) * MIntZ(It,Ie-1,IJ)
            END DO
            DO IJ = 1, NPrim
               MIntX(It,Ie,IJ) = MIntX(It,Ie,IJ) + ExpPHalf(IJ) * MIntX(It+1,Ie-1,IJ)
               MIntY(It,Ie,IJ) = MIntY(It,Ie,IJ) + ExpPHalf(IJ) * MIntY(It+1,Ie-1,IJ)
               MIntZ(It,Ie,IJ) = MIntZ(It,Ie,IJ) + ExpPHalf(IJ) * MIntZ(It+1,Ie-1,IJ)
            END DO
         END DO
      END DO

   END SUBROUTINE fmm_build_boxed_Mints

!-------------------------------------------------------------------------------
!
!   SUBROUTINE fmm_build_Mints(batch,NPrim,MaxAng,MaxMul)
!
!!     o Evaluate 1D multipole-moment integrals
!
!      IMPLICIT NONE
!
!      TYPE(fmm_prim_batch), INTENT(IN) :: batch(:)
!      INTEGER(INTK),        INTENT(IN) :: NPrim, MaxAng, MaxMul
!
!      INTEGER :: i, Ie, It, IJ
!
!      REAL(REALK) :: PXCX(NPrim), PYCY(NPrim), PZCZ(NPrim)
!      REAL(REALK) :: ExpPHalf(NPrim), ExpntP(NPrim)
!
!      PXCX(:) = batch(1:Nprim)%PC(1)
!      PYCY(:) = batch(1:Nprim)%PC(2)
!      PZCZ(:) = batch(1:Nprim)%PC(3)
!      ExpPHalf(:) = batch(1:Nprim)%ExpPHalf
!      ExpntP(:) = batch(1:Nprim)%ExpntP
!
!!     o Recurence relation
!!       M(t,e,NPrim)
!
!      MIntX(:,:,:) = Zero
!      MIntY(:,:,:) = Zero
!      MIntZ(:,:,:) = Zero
!
!      MIntX(0,0,1:NPrim) = SQRT(Pi/ExpntP(:))
!      MIntY(0,0,1:NPrim) = SQRT(Pi/ExpntP(:))
!      MIntZ(0,0,1:NPrim) = SQRT(Pi/ExpntP(:))
!
!      DO Ie = 1, MaxMul
!         DO It = 0, Ie
!            DO IJ = 1, NPrim
!               MIntX(It,Ie,IJ) = DBLE(It) * MIntX(It-1,Ie-1,IJ)
!               MIntY(It,Ie,IJ) = DBLE(It) * MIntY(It-1,Ie-1,IJ)
!               MIntZ(It,Ie,IJ) = DBLE(It) * MIntZ(It-1,Ie-1,IJ)
!            END DO
!            DO IJ = 1, NPrim
!               MIntX(It,Ie,IJ) = MIntX(It,Ie,IJ) + PXCX(IJ) * MIntX(It,Ie-1,IJ)
!               MIntY(It,Ie,IJ) = MIntY(It,Ie,IJ) + PYCY(IJ) * MIntY(It,Ie-1,IJ)
!               MIntZ(It,Ie,IJ) = MIntZ(It,Ie,IJ) + PZCZ(IJ) * MIntZ(It,Ie-1,IJ)
!            END DO
!            DO IJ = 1, NPrim
!               MIntX(It,Ie,IJ) = MIntX(It,Ie,IJ) + ExpPHalf(IJ) * MIntX(It+1,Ie-1,IJ)
!               MIntY(It,Ie,IJ) = MIntY(It,Ie,IJ) + ExpPHalf(IJ) * MIntY(It+1,Ie-1,IJ)
!               MIntZ(It,Ie,IJ) = MIntZ(It,Ie,IJ) + ExpPHalf(IJ) * MIntZ(It+1,Ie-1,IJ)
!            END DO
!         END DO
!      END DO
!
!   END SUBROUTINE fmm_build_Mints
!
!-------------------------------------------------------------------------------

   SUBROUTINE fmm_build_boxed_Mpole_ints(IAnglA,IAnglB,NPrim,MaxMul)

      IMPLICIT NONE

      INTEGER(INTK), INTENT(IN) :: IAnglA, IAnglB, NPrim, MaxMul

      INTEGER(INTK) :: It1, It2, Ie, It, IJ

!      MpoleX(0:IAnglA,0:IAnglB,0:MaxMul,1:NPrim) = Zero
!      MpoleY(0:IAnglA,0:IAnglB,0:MaxMul,1:NPrim) = Zero
!      MpoleZ(0:IAnglA,0:IAnglB,0:MaxMul,1:NPrim) = Zero
      MpoleX(:,:,:,:) = zero
      MpoleY(:,:,:,:) = zero
      MpoleZ(:,:,:,:) = zero

      DO It1 = 0, IAnglA ! s,p,d,f,,,
         DO It2 = 0, IAnglB
            DO Ie = 0, MaxMul
               DO It = 0, It1+It2
                  DO IJ = 1, NPrim
                     MpoleX(It1,It2,Ie,IJ) = MpoleX(It1,It2,Ie,IJ) +    &
                                             ECoefX(It1,It2,It,IJ) *    &
                                             MIntX(It,Ie,IJ)
                     MpoleY(It1,It2,Ie,IJ) = MpoleY(It1,It2,Ie,IJ) +    &
                                             ECoefY(It1,It2,It,IJ) *    &
                                             MIntY(It,Ie,IJ)
                     MpoleZ(It1,It2,Ie,IJ) = MpoleZ(It1,It2,Ie,IJ) +    & 
                                             ECoefZ(It1,It2,It,IJ) *    &
                                             MIntZ(It,Ie,IJ)
                  END DO
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE fmm_build_boxed_Mpole_ints

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_build_boxed_SphMpole(basis,batch,Ish,Jsh,NPrim,MaxMul)

      USE fmm_car_to_sph, ONLY: fmm_transform_car_to_sph

      IMPLICIT NONE
      TYPE(fmm_basis),      INTENT(IN) :: basis
      TYPE(fmm_prim_batch), INTENT(IN) :: batch(:)
      INTEGER(INTK),        INTENT(IN) :: Ish, Jsh, NPrim, MaxMul

      LOGICAL :: IEqJ
      INTEGER :: IAnglA, IAnglB
      INTEGER :: IL1, IL2, IL3
      INTEGER :: IeX, IeY, IeZ, Ie, It, MType, MCar, MSph
      INTEGER :: It1, Iu1, Iv1
      INTEGER :: It2, Iu2, Iv2
      INTEGER :: IJ, IBatch
      INTEGER :: IL2Temp
      INTEGER :: ndim
      REAL(REALK) :: Temp1, Temp2

      IAnglA = basis%KType(Ish)
      IAnglB = basis%KType(Jsh)

      ndim =  (basis%LtuvMax_Car(IAnglA) - basis%LtuvMin_Car(IAnglA) + 1)   &
            * (basis%LtuvMax_Car(IAnglB) - basis%LtuvMin_Car(IAnglB) + 1)
!FIXME do we need the next line??
      ndim = ndim * NPrim

      ALLOCATE(CarMpole(ndim, (MaxMul+1)*(MaxMul+2)/2, 0:MaxMul))
      IF (.NOT.ALLOCATED(SphMpole)) THEN
         ALLOCATE(SphMpole(ndim, (2*MaxMul+1), 0:MaxMul))
      END IF

      CarMpole(:,:,:) = zero

      IEqJ = (ISh == JSh)
      IBatch = 0
      DO IL1 = basis%LtuvMin_Car(IAnglA), basis%LtuvMax_Car(IAnglA)
         It1 = basis%Lt(IL1)
         Iu1 = basis%Lu(IL1)
         Iv1 = basis%Lv(IL1)

         IL2Temp = basis%LtuvMax_Car(IAnglB)
         IF (IEqJ) IL2Temp = IL1
         DO IL2 = basis%LtuvMin_Car(IAnglB), IL2Temp
            IBatch = IBatch + 1
            It2 = basis%Lt(IL2)
            Iu2 = basis%Lu(IL2)
            Iv2 = basis%Lv(IL2)

            DO MType = 0, MaxMul
               DO MCar = 1, (MType+1)*(MType+2)/2
                  IL3 = MCar + basis%LtuvMin_Car(MType) - 1
                  IeX = basis%Lt(IL3)
                  IeY = basis%Lu(IL3)
                  IeZ = basis%Lv(IL3)

                  Temp1 = Zero
                  DO IJ = 1, NPrim
                     Temp2 = MpoleX(It1,It2,IeX,IJ) *    &
                             MpoleY(Iu1,Iu2,IeY,IJ) *    &
                             MpoleZ(Iv1,Iv2,IeZ,IJ) *    &
                             batch(IJ)%CCoefAB
                     Temp1 = Temp1 + Temp2
                  END DO

                  CarMpole(IBatch,MCar,MType) = Temp1

               END DO
            END DO

         END DO
      END DO

      ! Transform batch of Cartesian integrals to spherical moments
      ndim = SIZE(CarMpole,1)
      CALL fmm_transform_car_to_sph(CarMpole,SphMpole,ndim,MaxMul)

      DEALLOCATE(CarMpole)

   END SUBROUTINE fmm_build_boxed_SphMpole

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_store_boxed_SphMpole(basis,batch,Ish,Jsh,MaxMul,empty_box)

      IMPLICIT NONE
      TYPE(fmm_basis),      INTENT(IN)  :: basis
      TYPE(fmm_prim_batch), INTENT(IN)  :: batch(:)
      INTEGER(INTK),        INTENT(IN)  :: ISh, JSh, MaxMul
      LOGICAL,              INTENT(OUT) :: empty_box

      REAL(REALK), PARAMETER :: MomScrn = 1d-15
      LOGICAL :: IEqJ
      INTEGER :: IAnglA, IAnglB
      INTEGER :: IL1, IL2, IJ
      INTEGER :: MType, MSph
      INTEGER :: IBatch
      INTEGER :: IL2Temp
      INTEGER :: Labelp, Labelq

      IEqJ = (ISh == JSh)
      IAnglA = basis%KType(Ish)
      IAnglB = basis%KType(Jsh)

      IBatch = 0
      Labelp = basis%KLoc_Car(ISh)
      DO IL1 = basis%LtuvMin_Car(IAnglA), basis%LtuvMax_Car(IAnglA)
         Labelp = Labelp + 1
         Labelq = basis%KLoc_Car(JSh)
         IL2Temp = basis%LtuvMax_Car(IAnglB)
         IF (IEqJ) IL2Temp = IL1
         DO IL2 = basis%LtuvMin_Car(IAnglB), IL2Temp
            IBatch = IBatch + 1
            Labelq = Labelq + 1

            DO MType = 0, MaxMul
               DO MSph = 1, (MType+MType+1)
                  ! Weak screening
                  IF (ABS(SphMpole(IBatch,MSph,MType)) < MomScrn) CYCLE
                  empty_box = .FALSE. 
                  WRITE(IO2) Labelp, Labelq, MType,     &
                             (MSph-MType-1), SphMpole(IBatch,MSph,MType)
               END DO
            END DO
            WRITE(IO2) -1, 0, 0,0, 0d0  ! marker for new (p,q)

         END DO
      END DO

      IF (ALLOCATED(SphMpole)) DEALLOCATE(SphMpole)

   END SUBROUTINE fmm_store_boxed_SphMpole

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_pack_boxed_mpoles(lmax,nBas,dens,nmoms_out)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: lmax, nBas
      REAL(REALK),   INTENT(IN)  :: dens(nBas,nBas)
      INTEGER(INTK), OPTIONAL, INTENT(OUT) :: nmoms_out

      REAL(REALK),       PARAMETER :: MomScrn = 1d-15
      CHARACTER(LEN=10), PARAMETER :: Name = 'multipoles'
      CHARACTER(LEN=255) :: FBuf
      REAL(REALK)   :: mpole(-lmax:lmax,0:lmax)
      REAL(REALK)   :: box(3), sph, T0, T1
      INTEGER(INTK) :: IOread, IOwrite1, IOwrite2
      INTEGER(INTK) :: i,j, l,m, nbox
      LOGICAL       :: empty_box

      CALL CPU_TIME(T0)
      IOread   = LUINTM
      IOwrite1 = LUINTM +1
      IOwrite2 = LUINTM +2

      ! File of AO-pair moments
      FBuf = TRIM(Name)//".fmm0"
      OPEN(UNIT=IOread, FILE=TRIM(FBuf), STATUS='OLD',ACTION='READ',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IOread)

      ! Open file of AO-pair moments summed over boxes
      FBuf = TRIM(Name)//".fmm1"
      OPEN(UNIT=IOwrite1, FILE=TRIM(FBuf), STATUS='REPLACE',    &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IOwrite1)

      ! Open file of box centres at which the FMM potential is computed
      FBuf = TRIM(Name)//".fmm2"
      OPEN(UNIT=IOwrite2, FILE=TRIM(FBuf), STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IOwrite2)

      WRITE(IO1,*) 'Writing file of boxed multipole moments...'
      nbox = 0
      box_loop: DO
         mpole(:,:) = zero
         shell_loop: DO
            READ(IOread) i,j, l,m, sph
            IF (i == -2 .AND. j == -2) EXIT box_loop
            IF (i == -1 .AND. j == -1) EXIT shell_loop
            IF (i < 0) CYCLE  ! marker for new (i,j)
            mpole(m,l) = mpole(m,l) + sph*dens(i,j)
         END DO shell_loop

         READ(IOread) empty_box, box(:)
         IF (.NOT. empty_box) THEN
            nbox = nbox +1
            ! Write summed moments to file
            DO l = 0, lmax
               DO m = -l, l
                  WRITE(IOwrite1) nbox,l,m,1,1, box(:),mpole(m,l)
               END DO
            END DO
            ! Write LHS grid point data
            WRITE(IOwrite2) 0,0,0,0,0, box(:), 1.0D0
         END IF
      END DO box_loop
      CLOSE(UNIT=IOread, STATUS='KEEP')

      ! Mark end of RHS moments file with negative L-index and close files
      WRITE(IOwrite1) 0,-1,0,0,0, 0d0,0d0,0d0, 0d0
      CLOSE(UNIT=IOwrite1, STATUS='KEEP')
      CLOSE(UNIT=IOwrite2, STATUS='KEEP')

      ! Write RHS boxed moments header file
      FBuf = TRIM(Name)//".fmm1header"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='REPLACE',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      WRITE(IO2) lmax, nBas, nbox
      CLOSE(UNIT=IO2, STATUS='KEEP')

      ! Write LHS potential grid points header file
      FBuf = TRIM(Name)//".fmm2header"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='REPLACE',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      WRITE(IO2) nbox
      CLOSE(UNIT=IO2, STATUS='KEEP')

      print *, 'Filled boxes, after packing =', nbox
      IF (PRESENT(nmoms_out) ) nmoms_out = nbox

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED in fmm_pack_boxed_mpoles', T1-T0, IO1)

   END SUBROUTINE fmm_pack_boxed_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_contract_multipole_potential(nBas,nbox,lmax,V,Jmat)

      USE fmm_qlm_utils, ONLY: fmm_renormalise_qlm

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)    :: nBas, nbox, lmax
      REAL(REALK),   INTENT(IN)    :: V(:,:)
      REAL(REALK),   INTENT(INOUT) :: Jmat(nBas,nBas)

      CHARACTER(LEN=10), PARAMETER :: Name = 'multipoles'
      CHARACTER(LEN=255) :: FBuf
      REAL(REALK)   :: mpole((1+lmax)**2,1)
!      REAL(REALK)   :: rho(nbox)
      INTEGER(INTK) :: n,i,j,ii,jj, l,m,lm, ibox
      REAL(REALK)   :: nel, sph, box(3)
      LOGICAL :: empty_box, screened_away

!print *, 'Potential:'
!do i = 1, size(V,2)
!   print *, i, V(1,i) 
!End do

      ! File of AO-pair moments
      FBuf = TRIM(Name)//".fmm0"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='OLD',ACTION='READ',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

      WRITE(IO1,*) 'Contracting file of raw moments...'
!      rho(:) = zero
      ibox = 1
      box_loop: DO
         mpole(:,1) = zero
         ii = 1; jj = 1 
         screened_away = .TRUE.
         shell_loop: DO
            READ(IO2) i,j, l,m, sph
            IF (i == -2 .AND. j == -2) EXIT box_loop
            IF (i == -1 .AND. j == -1) EXIT shell_loop
            IF (i<0 .AND. (.NOT. screened_away)) THEN
               ! Renormalise moments before contraction with potential
               CALL fmm_renormalise_qlm(lmax,mpole) 
               ! Contract potential for this (i,j) in this non-empty box
!               rho(ibox) = rho(ibox) + mpole(1,1)  ! need to include dens(i,j)
               Jmat(ii,jj) = Jmat(ii,jj) + DOT_PRODUCT(V(:,ibox),mpole(:,1))
               Jmat(jj,ii) = Jmat(ii,jj)
               mpole(:,1) = zero
               screened_away = .TRUE.
            ELSE IF (i > 0) THEN
               screened_away = .FALSE.
               lm = l*(l+1) +m +1
               mpole(lm,1) = mpole(lm,1) + sph
               ii = i; jj = j
            END IF
         END DO shell_loop

         READ(IO2) empty_box, box(:)
         IF (.NOT. empty_box) ibox = ibox+1

      END DO box_loop
      CLOSE(UNIT=IO2, STATUS='KEEP')
      IF (ibox-1 /= nbox) CALL fmm_quit('nbox /= ibox')

!nel = zero
!do i = 1, nbox
!   print *, i, rho(i)
!   nel = nel + rho(i)
!end do
!print *, "Total density:", nel

   END SUBROUTINE fmm_contract_multipole_potential

!-------------------------------------------------------------------------------

END MODULE fmm_boxed_multipole_ints
